library(tidyverse)
library(plotly)  
library(DT)
library(limma)  
library(clusterProfiler)  
library(org.Hs.eg.db)  
library(enrichplot)  
library(WGCNA)
library(sva)
library(VennDiagram)
library(ggvenn)

1 LIMA

1.1 Rejection

rejection = readRDS("/Users/PengOlivia/Desktop/PROMAD-Download-Rejection-2024-08-17.Rds")
rejection = rejection$GSE98320

column_names <- colnames(rejection)
outcome = column_names[-c(1, length(column_names))]

rejection_vector <- ifelse(grepl("AR", outcome), 1, 0)

rejection_exprs_mat = as.matrix(rejection[, -c(1, 857)])

design <- model.matrix(~ rejection_vector)
limma.fit <- lmFit(rejection_exprs_mat, design)
limma.fit <- eBayes(limma.fit)
rejection_DEresults <- topTable(limma.fit, number = Inf, sort.by = "p")

rejection_DEresults$Gene <- rejection$X

DErejection <- rejection_DEresults[,c("Gene", setdiff(names(rejection_DEresults),"Gene"))]
datatable(DErejection)

1.2 Fibrosis

fibrosis = readRDS("/Users/PengOlivia/Desktop/PROMAD-Download-Fibrosis-2024-08-17 new.Rds")
fibrosis = fibrosis$GSE76882

column_names <- colnames(fibrosis)
outcome = column_names[-c(1, length(column_names))]

fibrosis_vector <- ifelse(grepl("Fibrosis", outcome), 1, 0)

fibrosis_exprs_mat = as.matrix(fibrosis[, -c(1, 143)])

design <- model.matrix(~ fibrosis_vector)
limma.fit <- lmFit(fibrosis_exprs_mat, design)
limma.fit <- eBayes(limma.fit)
fibrosis_DEresults <- topTable(limma.fit, number = Inf, sort.by = "p")

fibrosis_DEresults$Gene <- fibrosis$X

DEfibrosis <- fibrosis_DEresults[,c("Gene", setdiff(names(fibrosis_DEresults),"Gene"))]
datatable(DEfibrosis)

2 GSCE

2.1 Enrichment plot for all rejection, fibrosis genes

#rejection
rejection_vector = DErejection$Gene
fibrosis_vector = DEfibrosis$Gene
biomarkers = intersect(rejection_vector, fibrosis_vector)

gsea_input = DErejection %>% 
  filter(Gene %in% biomarkers)

gene_list <- -log10(gsea_input$P.Value)
names(gene_list) <- biomarkers

gse <- gseGO(geneList=gene_list, 
             ont ="BP", 
             keyType = "SYMBOL", 
             OrgDb = org.Hs.eg.db)

datatable(gse@result)
path.id = gse@result$ID
gseaplot(gse, path.id[1])

#fibrosis 
gsea_input = DEfibrosis %>% 
  filter(Gene %in% biomarkers)

gene_list <- -log10(gsea_input$P.Value)
names(gene_list) <- biomarkers

gse <- gseGO(geneList=gene_list, 
             ont ="BP", 
             keyType = "SYMBOL", 
             OrgDb = org.Hs.eg.db)

datatable(gse@result)
path.id = gse@result$ID
gseaplot(gse, path.id[1])

2.2 Venn diagram for top 5 overlapping genes

#venn diagram

rejection_vector_top_500 = DErejection[1:500, ]$Gene
fibrosis_vector_top_500 = DEfibrosis[1:500, ]$Gene

sets <- list(Set1 = rejection_vector_top_500, Set2 = fibrosis_vector_top_500)

ggvenn(
  sets,
  show_elements = F,
  fill_color = c("red", "green", "blue")
)